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ABSTRACT 

Low-mass disks orbiting a massive body can support "slow" normal modes, in which 
the eigenfrequency is much less than the orbital frequency. Slow modes are lopsided, 
i.e., the azimuthal wavenumber m = 1. We investigate the properties of slow modes, 
using softened self-gravity as a simple model for collective effects in the disk. We employ 
both the WKB approximation and numerical solutions of the linear eigenvalue equation. 
We find that all slow modes are stable. Discrete slow modes can be divided into two 
types, which we label g-modes and p-modes. The g-modes involve long leading and long 
trailing waves, have properties determined by the self-gravity of the disk, and are only 
present in narrow rings or in disks where the precession rate is dominated by an external 
potential. In contrast, the properties of p-modes are determined by the interplay of self- 
gravity and other collective effects. P-modes involve both long and short waves, and in 
the WKB approximation appear in degenerate leading/trailing pairs. Disks support a 
finite number — sometimes zero — of discrete slow modes, and a continuum of singular 
modes. 

Subject headings: celestial mechanics, stellar dynamics — stars: formation — galaxies: 
nuclei 

1. Introduction 

Disks orbiting massive bodies are found in many astrophysical systems. These may be composed 
of gas (surrounding protostars, compact objects in binary star systems, or black holes in galactic 
nuclei), dust (a component of protoplanetary disks), solid bodies (planetary rings and planetesimal 
disks), or stars (galactic nuclei). Normally the disk mass in these systems is much less than the 
central mass M, so that the gravitational potential is nearly Keplerian. A special feature of the 
Keplerian potential is that eccentric orbits do not precess. This feature arises from a 1:1 resonance 
between the radial and azimuthal frequencies, and implies that the evolution of the eccentricity 
distribution in an isolated Keplerian disk is determined by the collective effects in the disk (self- 
gravity, pressure, collisions, viscosity, etc.), no matter how small these may be. Thus we expect 
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that Keplerian disks may be capable of supporting a set of normal modes with eigenfrequency u 
much less than the Keplerian orbital frequency il{r) = {GM /r'^Y^'^ , and properties determined by 
weak collective effects in the disk. In the linear limit these modes have azimuthal wavenumber 
m = 1. We shall call them "slow" modes. 

The goal of this paper is to investigate slow modes in Keplerian disks. In this Section we 
provide an introduction and orientation based on the WKB approximation. We shall find that in 
many cases slow modes have scales comparable to the size of the disk, so that the WKB analysis 
must be supplemented by numerical calculations of large-scale modes. We soften the self-gravity of 
the disk to mimic the effects of pressure in fluid disks and velocity dispersion in collisionless disks, 
and derive the integral eigenvalue equation that describes large-scale slow modes in §11. We discuss 
the numerical solution of this integral equation for several disk models in §111. Section IV contains 
a discussion of earlier related work, and §V contains conclusions. 



1.1. The precession rate 

We employ cylindrical coordinates (r, 0, z) with origin at the central mass M, and approximate the 
unperturbed disk as razor-thin, with surface density Sd(r) in the z = Q plane. We assume that 
the disk mass is small, M^jM ^ Y^^r^ jM ~ e 1. The unperturbed disk material orbits in the 
axisymmetric potential 

= + + (1) 

where ^aif) is the potential arising from the self-gravity of the disk and $e('') is the non-Keplerian 
potential from any external source, and both are assumed to be 0(e). 

For nearly circular orbits, the azimuthal and radial frequencies O > and k > are given by 
(e.g. Binney and Tremaine 1987) 

The line of apsides of an eccentric test-particle orbit subjected only to gravity precesses at a rate 



TX! = n — K 



1 +:^)($, + $e) + 0(6^), (3) 



2f2(r) \r dr dr"^ 
where from now on ri(r) is assigned the Keplerian value, (GM/r^)^/^. 
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1.2. The Kuzmin disk 

Several of our examples will be based on the Kuzmin disk, which has surface density and potential 
(e.g. Binney and Tremaine 1987) 

^'iW - 2;r(r2 + a2)3/2' " ~(^2 + „2)i/2' W 

where is the disk mass and a is the disk scale length. Using equation (3), the precession rate 
due to the Kuzmin disk is 

2^r)(r2 + a2)5/2' ^ ^ 

for M(i <ti M. The precession rate is negative at all radii, approaches zero as r^/^ for small 
radii and as r^''/^ at large radii, and has a single extremum at tq = {^Y^'^a = 0.6547a where 
w{ro) = -0.3257(Md/M)(GM/a3)V2. 



1.3. The WKB approximation 

We can write the linearized surface-density response of the disk as a superposition of terms of the 
form S(r, 0, t) = A{r) exp {i[J' k{r)dr + mcj) — wt]}, where m is the azimuthal wavenumber, k{r) 
is the radial wavenumber, and lo is the frequency. In the WKB or tight-winding approximation 
where |A;r| S> 1, the dispersion relation between k{r) and u depends only on the local properties of 
the disk, and can be derived analytically for both fluid and collisionless disks. 

For a barotropic fluid disk with sound speed c(r), the WKB dispersion relation reads (Safronov 
1960; Binney and Tremaine 1987) 

{u-mnf = K^ -2TTG^d\k\+k'^c^- (6) 

This relation implies that the disk is locally stable to axisymmetric (m = 0) disturbances if and 
only if 

For thin, low-mass disks we have c <C 0.r and <C M/r'^, so that in general the dispersion relation 
(6) can only be satisfied if \kr\ ^ 1, which justifies the use of the WKB approximation. However, 
the special case of m = 1 disturbances in a nearly Keplerian disk is different. Using equation (3) 
the dispersion relation can be rewritten in the form 

. TrGEJkl fcV 1^,.2 2^ 

Since we are interested in slow modes in nearly Keplerian disks, the terms of order w^, uP' can 
be dropped. At a given radius, the maximum value of uj is attained at wavenumber \k\ = ko = 
TrG^d/c^- 
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Let us neglect the effects of pressure for the moment, by dropping the term proportional to 
in equation (8). Thus the WKB dispersion relation for slow waves dominated by self-gravity 
becomes 

uj = w + . (9) 

Since oj/O = ©(E^r^/M) for isolated disks, in general the two terms on the right side of (8) are 
only comparable if \kr\ ~ 1; thus slow gravity-dominated waves must be large-scale (i.e. the WKB 
approximation is not satisfied) and have eigenfrequency uj = 0{w). The neglect of the pressure 
term is valid for such waves so long as GS^r ox M. ^ Q, where M.{r) = Vtr/c is the Mach 

number. 

Next consider a collisionless disk with a Schwarzschild distribution function having radial 
velocity dispersion Ci?(r). The WKB dispersion relation reads (Kalnajs 1965; Lin and Shu 1966; 
Binney and Tremaine 1987) 

(a; - mnf = - 27rGSd|A;|^ ( ^^^1^^ ^) , (lo) 

where ^(s, x) is the reduction factor. Explicit formulae for the reduction factor are given in the 
references above. The disk is locally stable to axisymmetric disturbances if and only if (Toomre 
1964) 

which closely resembles equation (7). For m = 1 disturbances in a nearly Keplerian disk, the 
dispersion relation (10) simplifies to 

u = w+ o ' ^ ^k\-^\. (12) 



in which 

^i^(x) = -e-^/i(x), (13) 
X 

and /i(x) is a modified Bessel function. At a given radius, the maximum value of is attained at 
wavenumbcr |A;| = /cq = 0.7643 ri/c^^. Since both terms on the right side of (12) scale as the disk 
surface density S^, for a fixed velocity-dispersion profile the shape of a mode is independent of the 
disk mass and the eigenfrequency is proportional to the disk mass. 

For the moment assume x ^ 1; then the reduction factor ^(x) = 1 and the dispersion relation 
(12) reduces to (9): we conclude, not surprisingly, that the WKB dispersion relation for slow 
gravity waves is the same for fluid and collisionless disks. The neglect of the reduction factor is 
valid for such waves so long as x 1 when |A;r| ~ 1, which in turn requires that the Mach number 
Mr = ^r/cR » 1. 

Many of the properties of waves in collisionless disks can be mimicked by softening the grav- 
itational potential, using —Gmjifi + 6^)^/^ instead of —Gm/d for the potential due to a mass m 
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at distance d. Here the parameter b is the softening length. This modification yields the WKB 
dispersion relation (Miller 1971; Erickson 1974; Toomre 1977) 

CO = w{r) + ^^^^1^1 exp(-|fc|6). (14) 

At a given radius, the maximum value of u is attained at wavenumber \k\ = = \/h. As in the 
case of a coUisionless disk, both terms on the right side scale as the disk surface density S^, so 
the eigcnfrequency is proportional to the disk mass. In the long-wavelength limit, \k\h <^ 1, the 
dispersion relation simplifies once again to equation (9). 

By analogy with the theory of WKB waves in galactic disks, we call waves trailing if A; > 
and leading if A; < 0. Long waves have \k\ < Aq and short waves have \k\ > k^. 

We conclude from these approximate arguments that a low-mass, thin disk (by which we mean 
E(^r^/M e ^ 1, and S> Q for a fluid disk, or S> 1 for a coUisionless disk) could support 
slow modes with the following properties: (i) the azimuthal wavenumber m = 1; (ii) the dominant 
collective effect is self-gravity; (iii) \kr\ = 0(1), that is, the wavelength is comparable to the disk 
radius, the properties of the mode are determined globally rather than locally, and the mode cannot 
be described accurately with the WKB approximation; (iv) the properties of the mode depend on 
the disk self-gravity but not on other collective effects; (v) the eigenfrequency is of order e$7. Many 
of these conclusions were derived already by Lee and Goodman (1999), in a paper focused mainly 
on nonlinear waves in nearly Keplerian disks. 

Despite conclusion (iii), the WKB approximation provides an invaluable guide to the properties 
of slow modes. It is useful to plot contours of constant u in the {k,r) plane. Figure 1 shows such 
contour plots for the Kuzmin disk discussed in §1.2 and four dispersion relations for slow waves: a 
disk with softened gravity (eq. 14), a fluid disk (eq. 8), a coUisionless disk (eq. 12), and a gravity- 
dominated disk (eq. 9). Local maxima in the contour plots occur at A; = ±A;o, and a local minimum 
at = 0. 

In these diagrams, a WKB wavepacket travels along a contour of constant to, at a radial speed 
given by the group velocity dto/dk (Toomre 1969). Long trailing and short leading waves propagate 
outwards, while long leading and short trailing waves propagate in. In the WKB approximation a 
mode corresponds to a closed contour in Figure 1 on which the total phase change § k{r)dr over 
one circuit is an integer multiple of 27r, including possible phase shifts at the turning points. The 
models in the Figure exhibit two types of closed contour: 

1. On contours such as C in the upper-left panel, an outward-traveling long trailing wave refracts 
into an inward-traveling short trailing wave, which in turn refracts back into a long wave. 
Contours such as B are similar, except that they involve inward-traveling long leading waves 
and outward-traveling short leading waves. WKB modes of types B and C occur in degenerate 
pairs. Modes of this kind are not present in the gravity-dominated disk (lower right panel), 
since this disk supports only long waves. In the Kuzmin disk these modes have a; > 0, and 
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Fig. 1. — Contours of constant frequency co for WKB slow waves in a Kuzmin disk (§1.2), as a 
function of wavenumbcr kr and radius r. The dispersion relations shown are for a disk with softened 
gravity (eq. 14), a fluid disk (eq. 8), a collisionless disk (eq. 12), and a gravity-dominated disk 
(eq. 9). The disk scalelength a = 1, the softening length b = 0.2r, the Mach number M = Cl/cr or 
Mr = ^r/cR is the same for the fluid and the collisionless disk and equal to 5 at all radii. Both 
the gravitational constant G and the central mass M are unity; the disk mass = 1 and u scales 
as so long as b, Mr or M^M'^ remain invariant. The extrema at kr = are minima. The 
contours in the right-hand panels are not uniformly spaced. 
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the modes with the smallest phase change around the contour have the highest frequency and 
wavemimbers localized near ifco- We call these p-modes, since they require the presence of 
pressure, velocity dispersion, or softening. 

2. In contrast, closed contours such as A involve only long waves, and involve reflection of an 
outward-traveling long trailing wave into an inward-traveling long leading wave, which then 
reflects back to long trailing. These modes are present in all four panels; their properties are 
similar in all panels since they depend only on the surface-density distribution in the disk 
and not other collective effects. In the Kuzmin disk these modes have a; < 0, and the modes 
with the fewest nodes have the most negative values of u. The WKB approximation is highly 
suspect for waves of this kind since \kr\ < 1, so numerical mode calculations are required. 
We call these g-modes, since self-gravity is the only collective effect involved. 

Figure 1 suggests that the WKB dispersion relation for collisionless Keplerian disks can be 
reproduced quite well by disks with softened gravity (Erickson 1974; Toomre 1977); and at least 
the qualitative features of the dispersion relation in fluid disks are similar as well. Therefore, in our 
numerical mode calculations we shall soften the self-gravity of the disk but neglect other collective 
effects. 



2. Linear perturbation theory 



In this Section we derive the linear eigenvalue equation for slow modes in disks with softened gravity. 
We write the surface density, velocity and potential of the perturbed disk as Sd(r) + Sa(r,t), 
Vd(r) + Va(r,t) = rO,(j) + Va(r,t) = rO0 + tta(r,t)r + Vair,t)4>, ^>d(r) + ^>a(r,t). The linearized 
Euler and continuity equations read 



dt 



+ (Vd • V)Va + (Va • V)vd 



dt 



0. 



(15) 



We write the perturbation variables {ua, Va,'^a, ^a) in the form Xa{r, (p, t) = X^{r) exp[i(m0— wt)]. 
Equations (15) become 



i{mn - a;)u^ - 20^^ 



dr ' 



+ - a;)C = -^<I.™, 

1 n 7777 

z(mn-a;)S- = --_(rS,n-) - — 



(16) 



The first two of these can be solved to yield 



I 

D 



d 2mQ 

{mil — CO)- — I 

dr r 
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v„ = 



1 

D 



— I (miz — u)) 

2n dr ' 



(17) 



where 



D = K^ - (mn-Lof. 



(18) 



2.1. Nearly Kepler ian disks 

We now specialize to slow perturbations in nearly Keplerian disks. Thus we assume that the disk 
potential and any external non-Keplerian potential satisfy l*!'^, ^e\/iGM/r) = 0(e) <ti 1, and that 
the azimuthal wavenumber m = 1. Furthermore we assume jwl/ri = 0(e), for the reasons given in 
§1.3. Using equation (3), the denominator D (eq. 18) becomes 



D = 2n{io - w) + 0{e^); (19) 



and equations (17) become 



^1 = i^(dK,^_K 

2{u -w)\dr r 

1 1 fd^l 2^l\ , , 

= 77 -7^ + — , 20 

" ^{uj-vj)\drr) ^ ' 

with fractional error 0(e). Thus to leading order in e, 

< = -2ivl, (21) 

which is simply the relation between the radial and azimuthal non-circular velocities for a nearly 
circular Keplerian orbit. The Keplerian eccentricity ex is related to the perturbed velocities by 

f GM\ ^^"^ 

v}^^ = —2iv\ = —i\ j exp(— iro). (22) 

We may use equation (21) to eliminate the radial velocity from the continuity equation (the 
last of eqs. 16); to leading order in e this now reads 

= ^-^ir^,vl) - ^ = 2^{^,vl) + (23) 
r dr r dr r 

The second of equations (20) can be rewritten as 

1 $^ 

Equations (23) and (24), together with Poisson's equation relating and describe the slow 
modes of Keplerian disks when self-gravity is the only collective effect. 
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2.2. Potential theory 

Poisson's equation can be written as 

$-(r)= / dr'r'PrnM^^ir'); (25) 
Jo 

the same relation, with m = 0, also holds between the potential and the surface density of 
the unperturbed disk. The kernel in equation (25) is 

Pmir,r') = -^6(-)(r</r>) + ^{d^, + S,,,.,). (26) 

Here r< = min (r,r'), r> = max(r, r'), and the second term is the indirect potential arising from 
the motion of the central body in response to the disk perturbation. The Laplace coefficient 

2 cos mQdO 



h/2i<^)-^l (i_2acos0 + a2)i/2 ^^^^ 



(Murray and Dermott 1999). In particular, 



bfha) = ^K{a), b'^ha) = ^[K{a) - E{a)], (28) 
' TT ' 7ra 



where K{a) and E{a) are complete elliptic integrals. We shall use the relations 

dE{a) _ E{a) - K{a) dK{a) _ E{a) _ K{a) 
da a ^ da a{l — a^) a 

The Laplace coefficients are singular at a = 1, 



(29) 



b\°) Jl -5) = ln5 + - ln8 + 0(51n(5); b^'^Ul -6) = ln5 + -(InS - 2) + 0{6lii6). (30) 

' TT TT ' TT TT 



We can soften the logarithmic singularity by replacing In 6 by ^ ln((^^ + P^) where /3 is the 
softening parameter; the actual softening length is then radius-dependent, b = f3r. If we adopt this 
prescription then we must soften the stronger singularities that arise from derivatives of hi 5 by 
differentiating \ \n{5'^ + /J^). Thus 

In5^iln(<52+A J -^4^, T^^WTW?- ^^'^ 
A full description of the softening procedure follows equations (45) and (53). 



2.3. The eigenvalue equation 



Combining equations (23), (24) and (25), we get 

r'dr' 



CO 



w{r)]vl{r) 



2^E,(/).iM + MMr01. (32) 
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This can be recast by writing 



0(r) 



1/2 



we then integrate by parts, using the fact that U{r) oc r to obtain 



f dr' 

ioy{r) =w{r)y{r) + I —K(r,r')y(r'), 



where the kernel 

K{r,r') = -2 



rSd(r)Ed(r')l 


1/2 ^ 







2 51n5 



2d\nr' 



Pi{r,r'). 



(33) 



(34) 



(35) 



The contribution of the indirect potential in equation (26) to K{r, r') vanishes, and the remaining 
component of Pi(r, r') is symmetric in r and r', so K{r,r') is also symmetric. 

Since the kernel is real and symmetric, the right side of equation (34) can be regarded as a 
linear, Hermitian operator acting on y{r'). Prom the properties of Hermitian operators we can con- 
clude (e.g., Courant and Hilbert 1953) that (i) the eigenvalues uj are real; thus all slow disturbances 
are stable; (ii) the eigenfunctions yn{r) associated with different eigenvalues Un are orthogonal with 
respect to the inner product {y, z) = J y*{r)z{r)dr/r, that is, (?/„, y^) = if a;„ 7^ w^. 

Since the kernel is real, the solutions y{r) or vl_{r) can be assumed to be real. In this case it 
proves useful to redefine the eccentricity ex (eq. 22) by the relation 



GM\ 



1/2 



exir); 



(36) 



so that ei^(r) is real but can have either sign. 

Using equations (26), (28) and (29) we can derive alternative forms for the kernel, 



K{r, r') 



ttG 



1/2 



= -G 



where a. = r^jr^ and 



C{a) = 2a-V2 



Sd(r)I]rf(rO 
rr'^{r)^{r') 

S,(r)S,(rO 
jr<^{r)^{r<) _ 

2(1 - + Q,4 



C(a), 



(l-a2)2 



-E{a) 



1 — a' 



K{a) 



(37) 



(38) 



The precession rate is given by equations (3) and (25) and can be written as 



w 



f dr' 



(39) 
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Here Wd and We are the precession rates from the self-gravity of the disk and from external sources 
(due to the potentials and $e respectively). The kernel is 



0(r)rV2 

where 



Dia), (40) 



D{a) = 2a^/^ 



l + a^ 



E{a) - :.K{a) 



.(l-a2)2 I- a 



(41) 



For the sake of simplicity we shall simply assume that any external precession mimics the 
precession rate due to the disk itself, that is, Weif) = f^di''^) where / is a constant. In this case 
the external precession can be eliminated if we replace the eigenfrequency uj by lo(1 + f) and replace 
the integral eigenvalue equation (34) by 

f dr' 

Loy{r) = Wd{r)y{r) + Xj -^K{r, r')y{r'), (42) 

where A = (1 + /)-^ 

The behavior of C{a) and D{a) as a is given by 

C{a) '-faV^ Dia) ^ (43) 

Both functions are singular at a = 1: if we write x = — log a > then 

^ 1 15, X 179 ^, , 

C{a) = — + — log- + — + 0(x), 



a;2 8 ° 8 48 

1 3, X 11 
^ + - log - + — 

x2 8 8 48 



D{a) = — + -log- + — + 0(x). (44) 



Figure 2 shows the behavior of C{a) and D(a); for clarity we have removed the dominant singularity 
by plotting c{a) = x'^C{a) and d{a) = x^d{a). 

The eigenvalue equation simplifies considerably in the case of narrow rings: we can replace 
C(a) and D(a) by the leading term in equations (44) and need not distinguish r from r' except 
where they occur in rapidly varying functions such as Sd(r) and w{r). Equations (33) and (39) 
become 

(7 f dr' 
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TTT^ 




0.0001 0.001 0.01 0.1 

1-a 

Fig. 2. — The functions C(a) and D{a) (eqs. 38 and 41). The dominant singularity has been 

removed by plotting c{a) = x'^C(a) and d{a) = x'^d(a), where x = — logo; (note that a < 1 so 
X > 0). The plot also contains dashed lines that almost coincide with the solid lines. These are the 
approximations used in the numerical integration scheme; the quantities plotted are x'^wk{x) (eq. 
48) and x'^wq{x) (eq. 53). 
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where fl can be treated as constant. The singular kernel can be softened according to equation 
(31): 

[u-^{rU{r) = -g/d/£,(rO.^(rO ^^^;j;y^^^^, , 
G r (r - r')2 - 52 

A trivial solution of these equations is a; = 0, 'Ua(r) =constant. 



2.4. Numerical methods 

Wc use similar numerical methods to evaluate the integrals in equations (39) and (42) and discuss 
both integrals together. We approximate the integrals using an A^-point grid that is uniform in 
u = logr. The main complication is the singularity of the integrands when r is close to r'. To 
handle this singularity, we write the integral in (42) as 

/ —K(r,r)y(r)= / dvwK(y-u) — (47) 

J r' J wk{v-u) 

where the weight function is chosen to mimic the singular behavior of the kernel (eq. 44) 



1 15. Ixl 179 

Is" 



^K{x)^-^ + -\og'-^ + —, (48) 



Before proceeding further we must soften the kernel to smooth the singularity. We do this by 
softening the weight function, that is, we replace equation (47) by 



j —K{r,r')y{r') = / dvwK{v — u) 



wk{v — u) 

Following the prescription in equation (31) the softened weight function is 



K{e-,e^)y{e^) ^^^^ 



X 



2 



/?2 15, (x2 + /32)i/2 179 



^^(^) ^ (^2+^2)2 + Y log 1 + (50) 

It should be stressed that this softened kernel does not correspond to any simple force law between 
two particles. Its key advantages are that it is numerically convenient when deriving a quadrature 
rule (see below); it has the correct behavior for softened gravity as re ^ 0; and the softened kernel is 
still symmetric in r and r' so the operator remains Hermitian. The softening not only simplifies the 
numerical evaluation of the integral but also can be used to mimic the effects of velocity dispersion 
in collisionless disks (cf. Fig. 1). 

Now if y(r') is smooth, the expression z{v) = K{e'^, e^)y{e^) Iwk{v — u) in equation (47) should 
also be smooth, so to evaluate the contribution to the integral from the interval [uk, Ufc+i] we apply 
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a four-point quadrature rule of the form 

dvK{e\ e-)y(e-) = J] WjZk+j-i, (51) 

where Zk = z(uk), and z_i = ^at+i = 0. The weights Wj are chosen so that the integration is exact 
if z{v) is a cubic polynomial, and can be evaluated analytically using the integrals v^WK{v—u)dv, 
n = 0, . . . ,3 (Press et al. 1992). 

Once the integral in equation (42) has been replaced by a sum, solving the integral equation 
reduces to finding the eigenvalues and eigenvectors of a symmetric N x N matrix, which can be 
done by standard methods (Press et al. 1992). 

A similar strategy is used to evaluate the precession rate: we write the integral in equation 
(39) as 

j l^rQir.r) = J dvwQiv - ^)— (52) 



where 



wq{v - ■ 

1 3 Ixl 11 
wq{x) = ^ + 8 log -§-^48^ ^^^^ 

the term ^x"^ is added to eliminate a zero of wq{x) that would otherwise occur at \x\ = 1.6712, 
but has no effect near the singularity at x = 0. We then soften wq{x) to wq{x) and recompute the 
weights Wj for the quadrature rule. 



3. Slow modes of disk models 

In our numerical calculations we typically use grids with N = 400-800. All the significant digits 
that we quote should be correct in the limit N ^ oo. 

Numerical mode calculations are more difficult with unsoftencd gravity than with softened 
gravity: the diagonal matrix elements become large when the softening is much less than the grid 
size, so the accuracy is degraded by roundoff error. Typically the minimum practical softening 
parameter was /3 = 10~^. 

Our numerical method always yields N eigenvectors and eigenvalues, but many of these are 
singular (i.e. the amplitude is much larger at a single grid point than elsewhere). These correspond 
to the singular (van Kampen) modes that occupy the continuous region of the eigenvalue spectrum. 



3.1. The Kuzmin disk 

We use units in which G = M = a = 1 (cf. §1.2). We usually employ a grid that covers the range 
— 7 < logr < 5. To account for precession due to the mass inside the innermost grid point, we 
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divide the disk into "live" and "frozen" components: the frozen component is a Kuzmin disk with 

log a = — 5 and the same central surface density as the original disk, while the live component is the 
difference between the original and frozen disks. This approximation is legitimate if the amplitude 
of the normal mode is negligible for logr < —5. We have checked that our results are not sensitive 
to the size of the frozen disk. 

Since softened self-gravity is the only collective effect, the shape of a mode is independent of 
the disk mass and the eigenfrequency is proportional to disk mass, so we are free to set = 1 as 
well, even though the calculations are only valid in real disks if ^ M. 

We classify the modes as g-modes and p-modes, as described in §1.3. 



3.1.1. g-modes 

We first examine Kuzmin disks with A = 0.1 and /3 ~ (actually /3 = 10"^), corresponding to the 
case in which 90% of the precession is due to an external source. The reasons for examining this 
somewhat artificial case are that (i) these disks do not support p-modes, since gravity is unsoftened; 
(ii) g-modes arc accurately described by the WKB approximation since A ^ 1. Thus we can both 
isolate the g-modes and compare their behavior to the analytic discussion in §1.3. 

This disk exhibits a rich spectrum of discrete slow modes, which can be sorted by the number 
of nodes; the modes with < 8 nodes are shown in Figure 3 along with their eigenfrequencies. As 
expected from the WKB analysis, the modes in Figure 3 all have a; < 0, and the modes with the 
fewest nodes have the most negative values of uj. 

The eigenfrequencies can be predicted using the WKB approximation. The dispersion relation 
(9) can be modified to include a factor A that accounts for precession from an external source (cf. 
eq. 42), and rewritten as 

\Hr)\= ^ [^-^(^)]- (54) 

The frequency contour diagram analogous to Figure 1 is shown in the left panel of Figure 4. The 
g-modes correspond to the contours centered on fcr = 0, r = tq = 0.6547; they have turning points 
at radii r± defined by w(r±) = co. The eigenfrequencies satisfy the quantum condition that the 
total phase change around the contour is an integer multiple of 27r, 



2 



j ^ \k{r)\dr = 2im, n = 1,2,3,... (55) 

J r- 



(there are additional phase shifts of ±^vr at the turning points, but these contributions cancel; see, 
for example, Shu et al. 1990). The right panel of Figure 4 shows the eigenfrequencies of the modes 
with 0, 1, 2, or 3 nodes, as a function of A, along with the WKB approximation from equations 
(54) and (55) for n = 1, . . . , 4. We see that the WKB approximation predicts the eigenfrequency 
fairly well for A < 0.3. At larger A, the WKB predictions are seriously in error. 
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Fig. 3. — Slow g-modes in the Kuzmin disk with A = 0.1 and /3 = 0. The modes are sorted by the 
number of nodes. The normahzation is chosen so that / ej^{r)dlogr = 1, where the eccentricity 
eft-(r) is defined in equation (36). The number in each panel is the eigenfrequency, in units where 
G = M = Md = a=l. The eigenfrequency scales as {Md/M){GM/a^y/'^. 
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Fig. 4. — (Left) Contours of constant frequency co for WKB slow waves in a Kuzmin disk with 
A = 0.1 and /? = 0, as a function of wavenumber kr and radius r. The contours are plotted at 
intervals of 0.05; the heavy contour is a; = 0. The g-modes correspond to the closed contours. 
(Right) The eigenfrequencies of g-modes with 0, 1, 2, or 3 nodes (bottom to top) for the Kuzmin 
disk, as a function of A (solid circles and lines). The dashed lines show the predictions of the WKB 
approximation (eqs. 54 and 55). There are no slow g-modes for A > 0.6. 
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As A increases, the frequency of each g-mode dechnes in absolute value, until the mode ter- 
minates when its frequency reaches zero. There arc no g-modcs in a Kuzmin disk with unsoftened 
gravity for A > 0.6. In particular, an isolated Kuzmin disk (A = 1) docs not support g-modes. 
The reason for this result can be understood qualitatively. As A increases (i.e. the importance of 
external precession declines), k{r) declines as A~^, (eq. 54) so that the integral on the left side 
of equation (55) also declines. Eventually this integral is less than 27r even for the largest closed 
contour in Figure 1, at which point the disk cannot support any g-modes. 



3.1.2. p-modes 

Next we examine p-modes in isolated Kuzmin disks with softened gravity. Here the WKB analysis 
is useful even for isolated (A = 1) disks so long as f3 <^ 1, since the characteristic wavenumber is 
ko = 1/b = l/(3r. The WKB analysis predicts that p-modes have to > 0; that the modes with the 
fewest nodes have the largest u; that the modes occur in degenerate pairs; and that the quantum 
condition is 

rr+ 

2 [k+{r)-k-{r)]dr = 2Tr{n-^), n = 1,2,3,... (56) 

J r_ 

Here k^ < ko < k-^ are the two solutions of the dispersion relation (14); r_ and r+ arc the turning 
points at which fc_ = A;_|_ for a given w-contour; and the term ^ on the right side arises because 
there is a phase shift of ^tt at each turning point (Shu et al. 1990). 

Figure 5 shows the predictions of this quantum condition for the Kuzmin disk, as a function of 
the softening parameter P = b/r. As /3 ^ the eigenfrequencies diverge as /3~^, so we have plotted 
/3a;. The limit point is 

55/4 

hm (3u = = 0.093575. (57) 

We have also plotted the ten largest eigenfrequencies from numerical calculations of the modes. 
These agree very well with the WKB approximation for small (3, and even for /3 = 0.3 the WKB 
estimates of the eigenfrequencies arc accurate within ~ 20%. Figure 6 shows the shape of the 
fundamental (n = 1) modes as a function of p. In contrast to the g-modes, the fundamental 
p-modes can have many nodes, because they are modulated at the high spatial frequency ko. 



3.2. Gaussian rings 



In this Section we examine slow modes in rings, by which we mean disks whose surface density 
approaches zero as r ^ 00 and as r ^ 0. As a model for rings we take the surface-density 
distribution 



Sd(r) = — exp 
r 



2s2 



(58) 
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Fig. 5. — Eigenfrequencies of slow p-modes in Kuzmin disks with softened gravity. The units are 

G = M = Md = a = 1. The filled eircles and solid curves show numerieal solutions for the ten 
highest eigenfrequencies. The dashed lines show the predictions from the WKB approximation (56) 
for n = 1, . . . , 5. In this approximation the modes appear in degenerate leading/trailing pairs. We 
have plotted (3uj instead of lo to expand the vertical scale (cf . eq. 57) . 
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Fig. 6. — Slow p-modes in the isolated Kuzmin disk with softened gravity. The two modes with 
the highest frequency are shown for disks with /3 = 0.05, 0.1, 0.2, and 0.4. Note that for /? ^ 1 the 
frequencies are nearly degenerate. The normalization is chosen so that / e|-(r)dlogr = 1, where 
the eccentricity exir) is defined in equation (36). 
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Fig. 7. — The apsidal precession rate zu for Gaussian rings with surface-density distribution de- 
scribed by equation (58). 
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where the constant K is related to the disk mass by Md = (27r)3/2irsexp(-is2). The precession 
rate in these rings is shown in Figure 7. Like the Kuzmin disks, Gaussian rings exhibit a single 
minimum in the precession rate, but the central minimum is surrounded by shoulders in which the 
precession rate is positive (prograde). 

3.2.1. Narrow rings 

We showed at the end of §2.3 that there is a trivial slow mode for narrow rings, with zero eigen- 
frequency and constant velocity perturbation. We can extend this analytic result with numerical 
solutions of the eigenvalue equation for narrow Gaussian rings. 

We find that slow modes for isolated rings with unsoftened gravity are present only if the ring 
is rather narrow, s < 0.07. Only the fundamental mode is non-singular; these are shown in Figure 
8. For s < 0.03 the shape of the fundamental mode is nearly flat, as expected from the trivial 
solution. At larger s the amplitude in the outer parts of the ring grows more and more sharply, 
until for s > 0.07 we were unable to find any non-singular modes. The frequency of the fundamental 
mode declined from lu = 1.1074 at s = 0.01 to a; = 0.6405 at s = 0.07, in units where the ring mass 

3.2.2. Broad rings 

We have examined a sequence of broader Gaussian rings, with softening parameter (3 = 0.2 and 
width s = 0.5, 1, 2. The three highest- frequency modes for each ring are shown in Figure 9. These 
are all p-modes. 

Like the Kuzmin disk, isolated broad Gaussian rings with unsoftened gravity do not support g- 
modes. We have examined a sequence of Gaussian rings with s = 1, near-zero softening {(3 = 10~^), 
and increasing A. The fundamental g-mode is well-defined for A = 0.5 (Figure 10), but disappears 
by A = 0.59. 

3.3. The Jacobs-Sellwood ring 

Jacobs and Sellwood (1999) have reported two-dimensional A^-body simulations of a low-mass 
annular disk orbiting a central mass. They were able to establish long-lived m = 1 normal modes 
that survived with little sign of decay for over 1500 orbital periods. The surface-density distribution 
in the disk was 

^a{r) = -^[l-^{l-rf]''\ i < r < §, (59) 
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Fig. 8. — The fundamental slow modes for narrow Gaussian rings (eq. 58). The modes are 
normalized so that J e^(r)cilogr/ Jdlogr = 1, and labeled by s, the dispersion in logr. There are 
no slow modes for s > 0.07. 




Fig. 9. — Slow p-modcs in the isolated Gaussian ring (eq. 58) with softened gravity, /3 = 0.2. The 
panels are labeled with the width parameter s and the eigenfrequency. The scale in logr is the 
same in all panels but the origin is different. 
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Fig. 10. — The fundamental g-mode in a Gaussian ring with width s = 1 and A = 0.5. Gravity is 
unsoftcned (/3 = 10"*^). The frequency of the mode is a; = —0.1527 and is proportional to the disk 
mass, which is here chosen to be = 1. 
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and zero otherwise. In contrast to our other models, this disk has sharp edges at r± = |, 5, where 
the surface density goes to zero as oc |r — r±|^/^. To ensure that our numerical methods handle 
edge effects properly, we extend the ring to /r+ and r_//, / = 1.35, with a very low surface-density 
envelope, = W'^Md. 

The precession rate in the unperturbed disk (59) is shown in the left panel of Figure 11, for 
several values of the softening parameter /3. The cusps at the inner and outer disk edge r± arise 
because of the sharp edges. 

The right panel of Figure 11 shows the frequencies of p- modes as a function of the softening 
parameter. The ten modes with the largest frequency are shown. The modes occur in degenerate 
pairs for /3 <IC 1 and the degeneracy is gradually lifted as /? increases. Each mode terminates at 
a maximum softening parameter, which is shown by a filled circle. For /? > 0.229 there are no 
non-singular modes. 

Figure 12 shows the shape of the highest frequency (fundamental) p-mode for several values of 
the softening parameter. Note that for the larger softening parameters, the amplitude of the mode 
does not go to zero at the disk edge. 

4. Earlier work 

Sridhar, Syer and Touma (1999) recognized the importance of slow modes in nearly Keplerian disks, 
and both formulated and solved numerically the linear eigenvalue problem. However, they employed 
the equations of Laplace-Lagrange secular perturbation theory (Murray and Dermott 1999), which 
is designed for planetary systems and so approximates the disk as a set of non-intersecting rings. 
This approach does not give the correct answers for a continuous disk, even as the number of 
rings becomes arbitrarily large, because it assumes that the eccentricities are small compared to 
the distance between rings. For example, all of the cigcnfrequencies are non-negative in Laplace- 
Lagrange theory (Sridhar et al. 1999), whereas we have found negative eigenfrequencies in both 
the WKB approximation and numerical solutions. 

Lee and Goodman (1999) have investigated stationary linear and nonlinear m = 1 density 
waves in nearly Keplerian disks. They have derived the nonlinear dispersion relation in the tight- 
winding approximation, and found self-similar logarithmic spirals with zero frequency in the special 
case where the disk surface density ^^(r) oc r"^/^. In the linear limit, their self-similar spirals 
exhibit properties reminiscent of the results of this paper: isolated disks do not support long waves, 
and zero-pressure disks in an external potential only support waves if the parameter A from equation 
(42) satisfies A < [r(i)/r(|)]^/144 = 0.53214, similar to our finding that there are no g-modes in 
the Kuzmin disk or s = 1 Gaussian ring if A > 0.5-0.6. 

Sridhar and Touma (1999) have examined the slow dynamics of eccentric orbits in nearly 
Keplerian potentials with a non-axisymmetric potential perturbation. They found two general 
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families of orbits, lenses and loops, and suggest that a self-consistent lopsided disk around a point 
mass could be assembled from aligned loop orbits. Our results extend this conclusion by showing 
that aligned loops generally can support self-consistent slow g-modes, but only if the precession 
is dominated by an external non-Keplerian potential or the disk is a narrow ring. In essence, 
the aligned loop orbits discussed by Sridhar and Touma work too well: the non-axisymmetric 
gravitational field that they produce is too strong for a self-consistent slow mode. The progeny of 
aligned loop orbits, which Sridhar and Touma (1999) call librating loop orbits, can support self- 
consistent p-modes in isolated disks, but large-scale p-modes require that the velocity dispersion is 
large. 

Statler (1999) has constructed sequences of periodic orbits in the potential of a point mass plus 
a lopsided disk. He argues that the precession of a self-consistent lopsided disk must be prograde 
(w > 0), and that the eccentricity ei^(r) must have a node. We find that neither of these results is 
generally true for g-modes: for example, the fundamental g-mode in the Kuzmin disk with A = 0.1 
(Fig. 3) has no nodes and a; < 0. Nevertheless, Statler's arguments may well describe the essential 
dynamics of some lopsided disks. 

Goldreich and Tremaine (1979) investigated slow modes in planetary rings. The thrust of their 
investigation was quite different, for several reasons: there was strong differential precession due 
to an external source (the planet's quadrupole moment), so that uniform precession was sustained 
only by a strong eccentricity gradient [de/dloga ^ 1); the ring was very narrow [Aa/a ~ 10~^); 
the ring was sharp-edged, so that the "resonant cavity" involved reflection at the disk edges rather 
than at turning points in the dispersion relation; and finally, later work (Chiang and Goldreich 
2000) shows that the mode properties in narrow, dense, particulate rings are strongly affected by 
collisions near the disk edge. In contrast, the slow modes examined here should be insensitive to 
the details of the disk edges since they are restricted to the interior of the disk (except for some of 
the modes in the strongly softened Jacobs-Sellwood ring. Figure 12). 



5. Conclusions 

Disks orbiting massive bodies support slow modes; these are linear m = 1 normal modes whose 
eigenfrequency is proportional to the strength of collective effects in the disk, rather than to its 
characteristic radial or azimuthal frequencies. The most important collective effects are likely to be 
the self-gravity of the disk, and the velocity dispersion in collisionless disks or the pressure in fluid 
disks. The first two of these and to a lesser extent the third — can be approximated by softened 
self-gravity (Fig. 1). Slow modes arc important because these are likely to be the only large-scale 
or long-wavelength modes, and hence can dominate the overall appearance of the disk. Moreover, 
slow modes are relatively immune to damping by viscosity. Landau damping, collisions, or other 
dissipative mechanisms. 

Slow modes in which the disk self-gravity is the dominant collective effect are described by 
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the same equations used in the classic Laplace-Lagrange secular perturbation theory for planetary 

systems, except that (i) the number of planets A'' goes to infinity and the mass per planet m oc N~^; 
(ii) the Laplace-Lagrange equations are valid in the limit where the planets' radial excursions are 
much less than their radial separation, whereas in continuous disks the radial separation of the 
mass elements is smaller than the radial excursions. The simplest way to modify Laplace-Lagrange 
theory to account for this difference is to soften gravity. 

We have derived the eigenvalue equation (42) that describes slow modes in a nearly Keplerian 
disk with softened gravity. We have solved this equation numerically to investigate the properties 
of modes in the Kuzmin disk (eq. 4), the Gaussian ring (eq. 58), the Jacobs-Sellwood ring (eq. 
59), and a variety of other disks not reported here. Many of our results can be interpreted in the 
framework of the WKB approximation (§1.3). Our conclusions can be summarized as follows: 

1. All slow modes have real frequency, and thus are stable. 

2. There are two main types of discrete slow mode: (i) g-modes involve long leading and long 
trailing waves. Their properties are insensitive to the softening parameter P once it is suf- 
ficiently small. In general these are retrograde (eigenfrequency co < 0). The fundamental 
g-mode has zero nodes, the lowest eigenfrequency, and a scale comparable to the disk radius, 
(ii) The properties of p-modes depend on the softening parameter; in particular their char- 
acteristic wavenumber is A; ~ In the WKB approximation that is valid for /? ^ 1, 
p-modes come in degenerate pairs, one composed of short and long trailing waves and the 
other of short and long leading waves. In general p-modes are prograde (w > 0), and the 
fundamental p-mode has the largest co. 

3. Narrow rings have a trivial slow g-mode, CKir) =constant. We have examined a wide variety 
of broad rings and disks, but have found none that support g-modes when they are isolated. 
In the presence of an external non-Keplerian potential that supplies a fraction 1 — A of the 
total precession, a variety of disks appear to support g-modes if and only if A < 0.5-0.6. 

4. The p-modes are modulated at the characteristic spatial frequency ko = (3/r. Thus long- 
wavelength p-modes occur only in disks with substantial softening (cf. Figures 6, 9, and 
12). 

5. In general disks support a finite number of discrete normal modes and a continuum of singular 
modes. In some cases, such as the Jacobs-Sellwood ring with (5 > 0.229, there are no discrete 
normal modes. 

Of course, these models have many shortcomings: they are linearized; they use softened gravity 
as a crude proxy for velocity dispersion or pressure; and the softening law (eq. 50) was chosen for 

mathematical convenience and does not correspond to a simple interparticlc force law. The models 
also do not address the important questions of how rapidly slow modes damp (although the N-body 
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simulations of Jacobs and Sellwood 1999 suggest they can be very long-lived), or how large the 
disk-mass ratio M^/M can be in disks that support discrete slow modes. To decide these issues 
there is no substitute for careful A^-body simulations; however, our results provide a guide for 
interpreting these simulations. This paper also does not address the question of how slow modes 
are excited and whether long-lived slow modes actually play a significant role in the structure and 
evolution of a variety of astrophysical disks. 

I thank Kathryn Johnston for a number of stimulating early discussions of this problem, and 
Jerry Sellwood and Alar Toomrc for thoughtfT.il comments. This research was motivated and guided 
by the simulations of eccentric disks carried out by Vince Jacobs and Jerry Sellwood, and I thank 
them for sharing their results freely. This work was supported in part by NSF grant AST-9900316 
and NASA Grant NAG5-7310. 
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Fig. 11. — (Left) The apsidal precession rate for the Jacobs-Sellwood ring (eq. 59) for softening 
parameter /? = 0.01,0.05,0.1,0.2. (Right) Slow p-modes in the isolated Jacobs-Sellwood ring with 
softened gravity. Wc show the ten modes with the highest frequency. The modes occur in degenerate 
pairs when /3 <C 1. Each mode terminates at a maximum softening parameter, shown by a filled 
circle. The units are chosen so that G = M = — 1, and a; cx M^. 
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Fig. 12. — The p-modc in the Jacobs-Senwood ring with the highest frequency, for /3 = 
0.02, 0.06, 0.12, 0.22. The dashed hnes show the amphtude of the mode in the a low surface-density 
envelope (Sd = 10~®Md) that surrounds the ring in our numerical model. 



